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Stochastic field equation for the canonical ensemble of a Bose gas 
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We present a novel norm preserving stochastic evolution equation for a Bose field. Ensemble 
averages are quantum expectation values in the canonical ensemble. This numerically very stable 
equation suppresses high-energy fluctuations exponentially, preventing cutoff problems to occur. 
We present 3D simulations for an ideal gas in various trapping potentials focussing on ground state 
occupation numbers and spatial correlation functions for a wide range of temperatures above and 
below the critical temperature. Although rigorously valid for non-interacting Bosons only, we argue 
that weakly interacting Bose gases may also be amenable to this approach, in the usual mean-field 
approximation. 
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Ultracold quantum gases in traps are currently being 
investigated to a hitherto unknown precision and in a 
variety of circumstances [U, 0, Q . The fascinating exper- 
imental possibilities to manipulate relevant parameters 
such as trap geometry, temperature, particle number, 
and even interaction strength show that these gases are 
ideal quantum systems to investigate and revisit funda- 
mental concepts of many particle and statistical physics. 

In these experiments, after cooling, traps contain 
roughly a fixed number of particles. Thus, from a physi- 
cal point of view, a canonical ensemble is to be preferred 
over a grand canonical one. For these finite systems, 
different predictions for fluctuations (even in the ther- 
modynamic limit [J]) call for a canonical description. 

We present a stochastic evolution equation for a c- 
number field ili{x) such that quantum statistical expec- 
tation values in the canonical state can be replaced by 
an ensemble mean over these stochastic fields. As ap- 
plications, we focus on densities, ground state occupa- 
tion numbers, and (spatial) correlation functions which 
have recently been measured in impressive experiments 

We stress two crucial properties of our novel equa- 
tion: first, the noise is spatially correlated, preventing 
cutoff problems to occur. Secondly, the equation is norm- 
preserving, reflecting the canonical nature of our ensem- 
ble. Both these properties ensure a very stable numer- 
ical solution such that full 3D problems may be tack- 
led. Moreover, our equation may be implemented in po- 
sition space such that arbitrary trapping potentials may 
be treated without any difhculty, and, eventually, inter- 
actions may be taken into account. 

While constructed for an ideal gas, we do strongly be- 
lieve that these positive features of our stochastic field 
equation (SFE) will also be valuable for the interacting 
case. There are a number of approaches that establish 
SFEs for the grand canonical state of an interacting Bose 
gas in a trap [H [IH, [H, [11 [11 . We will comment on 
these equations and their relation to our result towards 
the end of this work. 

We consider an ideal gas of N particles in a trap 
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with single-particle Hamiltonian H ^ + ^^i^) ~ 



X] £fekfc)(£fel ^-i^d eigenenergies e^. For the determination 
fc 

of A^-particle mean values {. . .) n =tr(. . . p^) we start 
with the canonical density operator 

PN = -^e-^/'^^P^ (1) 

in second quantization with the corresponding energy 
H = ^efed^afe, the canonical partition function Z^r, 

and the projector P/v = J2 \{^k}){{nk}\ onto the 

A^-particle subspace. As usual, the number states are 

|{"^fe}) = n(^l;)"'°/V«/c!|0), where nk is the occupation 

fe 

number of the k-th eigenstate |efc). 

With the notation iplx) — {xl-ip) the desired SFE 
takes the form of a nonlinear, norm preserving stochastic 
Schrodinger equation, here in Stratonovich calculus [l5l | 

^~-J{A + ^)H K^S^B] mt 

which is the central result of this paper. In equ. ([2]), 
A is a phenomenological damping parameter. Reflect- 
ing a fluctuation-dissipation-relation, it appears both, in 
the damping term and, as a square root, in the fluctua- 
tions. Crucially, we introduce an operator kT = ^h/^t_i 
that depends on the real temperature kT, obviously a 
reference to the Bose occupation number. The complex 
random field d^{x,t) = {x\d^{t)) represents white noise 
with correlations {d£^{t)\x){x'\d£^{t)) = 5{x — x')dt. Note, 
however, that the white noise is always acted upon by 

the operator V kT. For energies H <^ kT the operator 
kT acts simply as the multiplication with the thermal 
energy kT. However, for energies H ^ kT, the fluctu- 
ations are exponentially suppressed, which is a crucial 
feature of our novel equation to which we will come back 
towards the end of the paper. Another way of looking at 
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this is that proper quantum statistics leads to spatially 
correlated noise preventing the occurrence of arbitrarily 
high momenta [1^1 • 

We found equation ^ by starting with the Glauber- 
Sudarshan P-representation [17,] of the Gaussian expo- 
nential in ll]), 

= J -^Pir,nm)mi (3) 



Coherent states = \ipo)\ipi) ■ ■ ■ \'4'k) ■ ■ ■ are used 

for all modes and = YKe""/'''^ - 1) • 

k 




impNim) = M^)''eM-wm/Ni with (^i^) = 

S^fc IV'feP (see (3)) allows us to express (normally- 
ordered) quantum correlation functions as phase space 
integrals, for instance 



(4) 

The weight functions are given by W^Ar('0*,V') = 
e-'^'f'^'f'^ f (V'*,V') with the normalization con- 
stant Cn = / ^W^Ar(V'*,i/')- Note that second (or 
higher) order correlations have to be determined using 
Wn-2 (or lower index) in expression (H]), while the Cn 
remains. 

We stress that it is necessary to distinguish carefully 
between the norm Af — {tp\ip) of the stochastic field 
evolving according to equ. ([2]) - which remains con- 
stant for all times - and the particle number N. Apart 
from the solution of our (rigorous) SFE the exact 
stochastic simulation of the weight function Wn requires 
a distribution of values for J\f. It turns out, however, 
that for particle numbers much larger than one the norm 
distribution becomes narrow enough so that for all the 
temperatures and observables of interest in this pa per a 
simulation with a single norm Af is sufficient (see |16|). 
Still, there is a surprising subtlety: The distribution of 
the norm M depends on the absolute value of the ground 
state energy Eq. It is for Eq = only that we have to 
chose (Af) = N. The liberty to choose other values for 
Eq (and other norms Af, accordingly) can be used with 
benefit to achieve faster convergence in the numerical im- 
plementation (see [16]). 

Two further remarks are called for. First, being a 

highly nonlinear equation, we are not surprised to find 

that it is possible to replace the ensemble mean by a time 

t 

average W (x)ip(x'))N = hni 7 f ds ip* (x, s)il;(x' , s) 

over a single realization iplx^s). Secondly, we chose to 
propagate with the term (A + i) in our equ. ^ such 
that for A = the remaining complex unit describes 
"real" dynamics. In this way we can simulate the transi- 
tion from a non-equilibrium to an equilibrium state in a 
phenomenological manner - see recent experiments [l9| . 



FIG. 1: Momentum distribution of a Bose gas of 200 particles 
in a potential V{x, y, z) — x'^ + y'^ + 2^ during phase transi- 
tion. On the left hand side we show single realizations of our 
stochastic field equation ^ , on the right hand side long-time 
averages over 15000 time steps. Top: an almost "classical" 
distribution for T > T^; Middle: a peak develops for T ki T^; 
Bottom: most particles are condensed for T < Tc- 




FIG. 2: Ground state occupation as a function of tempera- 
ture for an ideal 3D Bose gas of 200 particles in a harmonic 
isotropic trap (plus signs) and for 100 particles in a box poten- 
tial (crosses). We compare with an analytical calculation for 
the harmonic (in a low temperature approximation) [23| (full 
line) and the box potential [2^ (dashed line). We also dis- 
play results for the thermodynamic limit (harmonic: dashed- 
dotted line, box: dotted line). 



We now turn to applications. While a gas in a box is 
best treated in momentum space, a general implementa- 
tion of equ. ^ in position space is advantageous, since it 
can be adjusted easily to any trapping potential (and in a 
next step mean field atomic interactions may be included 
- see later). However, in position space, the generation of 

the correlated noise {x\ V kT\d^) is cumbersome. For the 
simulations presented here we use a Wigner-Weyl repre- 
sentation [l3| of the operator kT and consider only terms 
of lowest order in h. As the examples below show, this 
approximation is legitimate (for more details, see [iq). 
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The functioning of our equation is visualized in Fig. 
[TJ We simulate a 3D Bose gas containing 200 particles, 
trapped in a quartic potential V{x, y, z) — x'^+y'^+z'^ and 
determine the momentum densities n{px,Py) (the third 
momentum is integrated over). Three different tem- 
peratures are chosen: above (top), at (middle), and be- 
low (bottom) the critical temperature for Bose-Einstein 
condensation. While on the left hand side we display a 
single realization of eq. ([2]) after a certain propagation 
time, the right hand side shows time averages over 15000 
time steps. Obviously, on average we obtain the typical 
pictures for the transition to a Bose-Einstein condensate. 

The true quality of our equation is to be verified by 
calculating various characteristic quantities. As shown 
in Fig. [1] the numerical code in position representation 
allows us to treat the Bose gas in any trapping poten- 
tial. In order to make contact to previous results for 
the canonical ensemble, however, we restrict ourselves in 
the following to a 3D box and a 3D harmonic oscillator 
potential. More general considerations will be published 
elsewhere fl6|. 

First we show the ground state occupation in Fig. [51 
The results for the harmonic oscillator (with N = 200 
particles, plus signs) arc obtained by propagating cqu. 
(I2|) on a position grid. No use is made of the known 
spectrum and eigenfunctions. We compare with an ana- 
lytical approximation (full line) [20| . which is known to 
be in good agreement with exact results, and the thermo- 
dynamic limit (dashed-dotted line) . Next we show results 
for a 3D box {N — 100 particles, crosses) obtained from 
our SFE (here computed in momentum space) compared 
with an (approximate) result based on a path integral ap- 
proach |2l| (dashed line) and find very good agreement. 
Temperature is scaled to the critical temperature of the 
thermodynamic limit [l|, Q . Note that finite size effects 
are very significant for the box as seen when comparing 
with the thermodynamic limit (dotted line). 

Next we determine spatial correlation functions of first 
Gi{x,x') = {tp'^{x)^{x'))N and second order G2{x,x') — 
{'ip'^ {x)%l:i'^ {x')4'{x)4'{x')) N above and below the critical 
temperature. Gi can be measured through interference 
experiments (see Q). G2 and related quantities have 
been measured more recently in impressive experiments 
[1, 0, H, Q ■ For first order correlations, the canonical re- 
sults are close to the grand canonical values. For second 
order correlations and temperatures below the critical 
temperature, however, large deviations appear and the 
results of the canonical ensemble may be obtained ap- 
proximately with the help of condensate corrections to 
the grand canonical calculation [22]. In Fig. [31 we show 
Gi(x,0) (top) and 62(2;, 0) (bottom) obtained from the 
SFE ([2]); here for a 3D Bose gas of 1000 particles in a box 
with periodic boundary conditions (left hand side) and 
for a gas of 200 particles in an isotropic harmonic oscil- 
lator (right hand side). Our values (plus signs, crosses, 
stars) are compared with a direct grand canonical cal- 
culation ("gc" in Fig. [3]) and the corrected description 
("gc-fc" in Fig. [3]) based on the theory of f22l|. These 
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FIG. 3: First and second order spatial correlations Gi(a;,0) 
(top) and G2(a;,0) (bottom) for a Bose gas in a 3D box po- 
tential with periodic boundary conditions (left) and a 3D har- 
monic oscillator potential (right) for different temperatures. 
Results from the SFE (plus signs, crosses, stars) are com- 
pared to a direct calculation based on the theory of 22]: a 
grand-canonical description ("gc", solid, dashed and dashed- 
dotted lines) fails for G2, while a corrected theory ("gc-(-c" 
dashed-dotted and dotted lines) shows good agreement. 



corrected second order correlation functions are in good 
agreement with our exact canonical calculations. 

We see the main achievement of this paper in the 
fact that we established a numerically very robust (ex- 
act) SFE for the canonical state of an ideal Bose gas, 
amenable to an efficient numerical solution in full 3D 
for arbitrary trapping potentials. Still, in the light 
of the wealth of activities involving ultracold quantum 
gases, it is certainly of great importance to also in- 
vestigate the interacting case. Let us therefore relate 
equ. ([2]) to previous stochastic equations constructed for 
the grand canonical ensemble of an interacting Bose gas 
[13, In!, [H, [H, [11. A good survey over these differ- 
ent findings, their relations and their limitations is given 
in Il3i. Most notably, in several of these approaches 
[iSliE HI, [11, termed "classical field methods" in [l3| 
(and see also [23|), due to ultraviolet problems, lowly oc- 
cupied states must be cut off fll] or treated in a different 
formalism [13. [T3j. 

Let us now turn to our SFE omitting the non- 
linear terms and substituting both H H — H — fiN 
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and kT ^kT = H/{e"/^'^-l) {n the chemical potential 
and N the number operator), equ. ^ reduces to 

d\^) = -(A + i)H\^) + \] 2AfcT \dCi (5) 

of which one can easily show that it indeed provides 
proper grand canonical ensemble averages [l^l- We has- 
ten to stress that the canonical equation ^ is not merely 
a normalized version of equ. ([5]). As discussed before, in 

(l5|) the operator kT incorporates a natural high energy 
noise cutoff induced by proper quantum statistics. No 
further care is required. 

Describing a grand canonical ensemble, equ. fS]) is 
the link to establish a connection to the "classical field 
methods" mentioned above for an interacting gas: First, 
the operator kT appears as the simple temperature kT in 
those approaches, requiring a cutoff. More importantly, 
the interaction may be taken into account by a mean field 
contribution H ~> H + g\'ip{x)\'^, where g = 4TTalT?/m is 
the interaction parameter and a the s-wave scattering 
length. 

After these considerations it appears more than tempt- 
ing to use equation ([2]) for a canonical ensemble even 
in the case of a weakly interacting gas, with V{x) — > 
V{x) + g\tp{x)\'^. As argued, the resulting equation is 
free from ultraviolet problems and coincides (in the grand 
canonical case) with previous ( "classical" ) equations (in- 



cluding interaction and a high-energy cutoff). Moreover, 
it reduces to the (imaginary-time) Gross-Pitaevskii equa- 
tion for kT — > 0, describing a pure condensate (with the 
given particle number N). 

Let us briefly summarize our result: we present a 
norm-preserving stochastic field equation for the canoni- 
cal state of a Bose gas, describing experiments with a fi- 
nite number of atoms in an arbitrary trap. Being driven 
by spatially correlated noise, cut-off issues do not ap- 
pear; the equation is numerically very stable. We stress 
that it is valid for arbitrary temperatures: for T ^ Tj, it 
provides a wave description of a classical gas of massive 
particles. This is very much in the spirit of the way we 
think of light emerging from a light bulb as being com- 
posed of incoherent wave trains. We are able to deter- 
mine important quantities like spatial correlation func- 
tions and occupation numbers as a function of temper- 
ature in arbitrary traps. Finally, we relate our equation 
to stochastic Gross-Pitaevskii equations that exist for in- 
teracting gases in the grand canonical ensemble, raising 
expectations that the new equation should be applicable 
to weakly interacting Bose gases as well. 
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port by the International Max Planck Research School 
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Dresden. 
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